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Proposed methods for measuring the electric dipole moment (EDM) of the proton use an intense, 
polarized proton beam stored in an all-electric storage ring “trap”. At the “magic” kinetic energy 
of 232.792 MeV, proton spins are “frozen”, for example always parallel to the instantaneous particle 
momentum. Energy deviation from the magic value causes in-plane precession of the spin relative 
to the momentum. Any non-zero EDM value will cause out-of-plane precession—measuring this 
precession is the basis for the EDM determination. A proposed implementation of this measurement 
shows that a proton EDM value of 10“^® e-cm or greater will produce a statistically-significant, 
measurable precession after multiply-repeated runs, assuming small beam depolarization during 
1000 second runspp, with high enough precision to test models of the early universe developed to 
account for the present day particle/anti-particle population imbalance. 

This paper describes an accelerator simulation code, ETEAPOT, a new component of the Unified 
Accelerator Libraries (UAL), to be used for long term tracking of particle orbits and spins in 
electric bend accelerators, in order to simulate EDM storage ring experiments. Though qualitatively 
much like magnetic rings, the non-constant particle velocity in electric rings give them significantly 
different properties, especially in weak focusing rings. Like the earlier code TEAPOT (for magnetic 
ring simulation) this code performs exact tracking in an idealized (approximate) lattice rather than 
the more conventional approach, which is approximate tracking in a more nearly exact lattice. The 
BMT equation describing the evolution of spin vectors through idealized bend elements is also 
solved exactly—original to this paper. Furthermore the idealization permits the code to be exactly 
symplectic (with no artificial “symplectification”). Any residual spurious damping or anti-damping 
is sufficiently small to permit reliable tracking for the long times, such as the 1000 seconds assumed 
in estimating the achievable EDM precision. 

This paper documents in detail the theoretical formulation implemented in ETEAPOT. The ac¬ 
companying paper[2] “EDM planning using ETEAPOT with a resurrected AGS Electron Analogue 
ring” describes the practical application of the ETEAPOT code in the UAL environment to “res¬ 
urrect”, or reverse-engineer, the “AGS-Analogue” all-electric ring built at Brookhaven National 
Laboratory in 1954. Of the (very few) all-electric rings ever commissioned, the AGS-Analogue ring 
is the only relativistic one and is the closest to what is needed for measuring proton (or, even more 
so, electron) EDM’s. That paper also describes preliminary lattice studies for the planned proton 
EDM storage rings as well as testing the code against analytic calculations. 

PACS numbers: 14.20.Dh, 29.20.Ba, 29.20.db, 29.90.-|-r, 42.25.Ja 


I. INTRODUCTION 

a. Orbit and spin simulation code needed for 
electric storage rings. The U.S. particle physics com¬ 
munity has recently updated its vision of the future and 
strategy for the next decade in a Particle Physics Project 
Prioritization Panel (P5) Report. One of the physics 
goals endorsed by P5 is measuring the EDM of funda¬ 
mental particles (in particular proton, deuteron, neutron 
and electron). 

Since Standard Model EDM predictions are much 
smaller than current experimental sensitivities, detec¬ 
tion of any particle’s non-zero EDM would signal dis¬ 
covery of New Physics. If of sufficient strength, such 
a source could provide an explanation for the observed 
matter/antimatter asymmetry of our universe. A pro¬ 
ton EDM collaboration^ has proposed a storage ring 
proton EDM measurement at the unprecedented level 
of 10“^®e- cm, an advance by nearly 5 orders of mag¬ 


nitude beyond the current indirect bound obtained using 
Hg atoms. 

This paper is limited to the theoretical orbit and spin 
dynamical formulation within ETEAPOT, which is a 
newly developed code within the Universal Accelerator 
Libraries (UAL) simulation environment. The accompa¬ 
nying paper “Using ETEAPOT to resurrect the AGS- 
Analogue ring for EDM planning” describes the practi¬ 
cal application of the ETEAPOT code with emphasis on 
details of simulation requirements for the EDM measure¬ 
ment. 

b. Complications imposed by electric bending. 

The fundamental complication of an electric ring, as con¬ 
trasted to a magnetic ring, is the non-constancy of par¬ 
ticle speed. A fast/slow separation into betatron and 
synchrotron amplitudes has become fundamental to the 
conventional (Courant-Snyder) magnetic ring formalism. 
But, in an electric lattice the mechanical energy (as quan¬ 
tified by the relativistic factor 7 = 1 /-^1 — ( 3 ^) varies on 



2 


the same time scale as the transverse x and y ampli¬ 
tudes. On the other hand, changing only in RF cavities, 
the total energy E = 'ymc^ + eV{r), which includes also 
the potential energy eV{r), changes on a slow time scale, 
which makes a similar fast/slow separation valid. 

In the magnetic formalism, since it is only in RF cav¬ 
ities that the mechanical energy varies, the 7 factor is 
invariant in the rest of the ring. Furthermore one is ac¬ 
customed to treating energy as constant for times short 
compared to the synchrotron period. To the extent the 
betatron parameters are independent of particle energy, 
the betatron and synchrotron motions can then be su¬ 
perimposed. 

To most closely mimic this treatment in an electric 
ring, and to continue to regard 7 as the fundamental 
“energy-like” parameter, requires us to evaluate 7 only 
in regions of zero electric potential, which is to say, not 
in RF cavities, and not in electric bending elements—in 
other words, only in field free drift regions. This leads 
to a curious, but entirely satisfactory, representation in 
which the particles spend most of their time inside bend 
elements where 7 is varying, and little time in short drift 
regions where 7 is constant. The reason this is fully sat¬ 
isfactory is that the drift regions are closely spaced, and 
more or less uniformly spaced around the ring. Know¬ 
ing the lattice functions exactly at those points is oper¬ 
ationally equivalent to knowing them everywhere. With 
this interpretation, one can again rely on the approxi¬ 
mate representation of the motion as a superposition of 
fast betatron and slow synchrotron motions. 


II. PARTICLE TRACKING PARADIGMS 

The conventional formulation of accelerator physics is 
based on a paraxial approximation in which all orbit an¬ 
gles are small relative to a central design orbit. Not only 
is this approximation quite good for small rings, it be¬ 
comes progressively better as ring radii increase. The 
most important formulas obtained in this approach are 
based on linearization of the paraxial orbit equations. For 
sufficiently small amplitudes orbit evolution can be repre¬ 
sented by “transfer matrices” multiplying initial displace¬ 
ments, represented in 6 D phase space by six component 
displacement vectors. By introducing nonlinear “trans¬ 
fer maps” this approach can be extended to larger am¬ 
plitudes by representing each output variable as a (trun¬ 
cated) power series of terms each of which is a product of 
initial components, or their squares, cubes, etc. This ap¬ 
proach can be referred to as approximate tracking in an 
“exact” lattice where “exact” means, for example, that 
magnetic field profiles, fringe fields etc., once represented 
faithfully by accurate power series, can be subsumed into 
the previously described truncated power series orbit rep¬ 
resentations. 

The TEAPOT tracking paradigm[3] is very different. 
It is a refinement of the “kick code” paradigm, under de¬ 
velopment when TEAPOT was first introduced, about 


thirty years ago. Though still based on the paraxial ap¬ 
proximation, elements are sliced into sufficiently short 
segments that they can be represented by delta function 
“kicks”, or kinks, where orbit slopes change discontinu- 
ously, but displacements are continuous. The main virtue 
of this approach is that it is symplectic; i.e. particle 
beams respect Liouville’s theorem. 

TEAPOT took the more extreme approach of jettison¬ 
ing the paraxial approximation altogether, and insisting 
that orbits be constructed only from exactly determined, 
analytic orbits, interupted where appropriate, by kicks. 
This can be referred to as exact tracking in an approxi¬ 
mate lattice. This approach is restricted by the fact that 
exact orbits are known only in very special cases; for ex¬ 
ample in uniform magnetic fields. However, by introduc¬ 
ing thin kick elements it is possible to represent actual 
elements and actual rings to good accuracy. These added 
kicks may be artificial, meaning they model field devia¬ 
tions from the idealized field, or they may represent the 
total fields of elements physically present in the lattice. 
All this is explained in detail in the original Schachinger 
and Talman paper[3], which remains valid, and in use, 
to this day. To restrict the length of the present paper, 
only features newly introduced in ETEAPOT, and not 
covered in the original paper, are described. 

One new feature is the requirement to use electric 
rather than magnetic bending elements, which brings in 
the complication already mentioned, that, unlike in mag¬ 
nets, particle speeds are not constant in electric elements. 
The other new feature is the need to also track the parti¬ 
cle spins. To respect the TEAPOT paradigm this further 
restricts the list of allowable lattice elements to elements 
for which the BMT equation, an abbreviation for Michel, 
Bargmann, and Telegdi[l], is exactly solvable. Fortu¬ 
nately it is possible, even straightforward, to meet both 
of these new restrictions. 

The ETEAPOT approach has been developed to ad¬ 
dress new requirements of the proton EDM project. First 
is the transition from mgnetic to electric bending. Mag¬ 
netic tracking codes like TEAPOT and MAD implicitly 
take advantage of the constant speed of particles in mag¬ 
netic fields. This seriously complicates the porting of 
simulation algorithms to electric rings in which the par¬ 
ticle speed depends on the local electric potential en¬ 
ergy. Thick bending elements are again to be replaced 
by integrable force fields. As it happens. Inverse square 
law is the unique force field for which the 3D relativistic 
equations of motion can be solved exactly. The second 
major ETEAPOT extension has been the introduction 
and tracking of particle spin vectors (which provides the 
EDM signature of the EDM experiment). The differen¬ 
tial equation governing spin motion is the BMT equation. 
In general 3D motion the BMT equations are computa¬ 
tionally challenging. But in our idealized elements the 
BMT equations can be solved exactly, as required by the 
fundamental TEAPOT/ETEAPOT design principle. 

By symmetry, every orbit in a radial force field lies 
in a plane, referred to here as “the bend plane”. This 
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plane is always close to but, in general, not quite identi¬ 
cal to the horizontal symmetry of the ring. In spherical 
(r, 9) coordinates the differential equation governing r(9) 
for a relativistic particle orbit in an inverse square law 
force field can be solved exactly [ 5 ]. These are like clas¬ 
sical planetary orbits or the orbits in a hydrogen atom 
treated classically. Even relativistically, every orbit is a 
processing ellipse-like (rosette) figure, familiar to Newton 
in classical mechanics, and for relativistic motion to Ein¬ 
stein, in his general relativistic calculation of the advance 
of the perihelion of Mercury. 


III. RELATIVISTIC KINEMATICS IN 
ELECTRIC POTENTIAL V{r). 

In the horizontal y=0 bend plane, a radial electric field 
with index m power law dependence on radius r is 

E(r,0) = -Eo^f, (2) 

and the electric potential V{r), adjusted to vanish at 
r = ro, is 


Using (r, 0) polar coordinates, it will be shown shortly 
that orbit evolution in an inverse square law bend plane 
is described exactly by 


r{0) 


A 

1 -I- e cos k{ 6 — 0o) ’ 


and X = r — ro 


( 1 ) 


where r is the radial distance from the bend center to 
the particle position. Then x is the (paraxial) radial dis¬ 
placement from the design central orbit, which is a cir¬ 
cle of radius rg. This formula is as simple as it is only 
because the motion is two dimensional. The potential 
energy depends only on r(9), permitting the speed, and 
hence the momentum components to be expressed ana¬ 
lytically as functions of 9. In correlating this equation 
with conventional paraxial treatment, one notes that the 
independent variable 9 here differs from the independent 
longitudinal variable s, which is the path length along the 
design orbit; but s and 9 advance proportionally in bend 
elements, with the constant of proportionality being the 
design radius rg. 

Each particle being tracked has its own private orbit 
plane parameters in the orbit equation—though they are 
all very nearly the same for realistic beams. The pa¬ 
rameters in Eq. Q are interpretable (approximately) as 
follows: A is “average” radius, e is “eccentricity”, 9o is the 
angular deviation from perihelion, and k (interpretable 
as the “tune” in accelerator jargon) establishes the beta¬ 
tron phase advance per unit angular advance. 

For particle tracking this formulation is simpler than 
the 3D description of orbits in a uniform magnetic field. 
One price to be paid for this simplicity is the need to 
transform each particle into, and later, out of, its private 
2D orbit plane (with its own basis vectors) as it enters, 
and later, leaves a bending element. These transforma¬ 
tions, not shown here, are elementary near-identity rigid 
rotations (including spin). 

A more serious price for this exact, “global” coordi¬ 
nate, approach is purely numerical. Though Eq. 0 is an¬ 
alytically exact it is numerically treacherous—the parax¬ 
ial quantity x is a millimeter scale length compared to 
r which is very nearly equal to the rg which is a large 
length, such as 40 m. The special numerical treatment 
needed to handle this issue is discussed below. 


U(r) 


Epro 

m 


( 3 ) 


The “cleanest” case, shown in Figure [l] has m=l and is 
known as the Kepler or the Coulomb electric field of a 
point source. One way our case is more general than Ke¬ 
pler’s is that our treatment has to be exactly relativistic. 
This m=l case can be referred to as “spherical” since the 
equipotentials are spheres centered on a point charge at 
the center, and r is the radius in (r, 9, (p) spherical coor¬ 
dinates. For m=l the Kepler problem can be solved in 
closed form with the same generality as in the relativistic 
as in the nonrelativistic case. However the orbits are no 
longer exactly elliptical, nor restricted to a single plane. 

The m=0 field is referred to as “cylindrical” and r is 
the radius in (r, p, y) cylindrical coordinates. This field 
is the experimentally easiest to produce, using cylindical 
electrodes having the same central axis. It has been tac¬ 
itly assumed that the EDM storage ring bends will be 
constructed this way, even if the optimal electrode shape 
would be somewhat different. As already explained, 
(weak) thin quadrupoles can be used to tailor the fields 
to overcome this impediment. Curved-planar electrodes 
produce the electric field shown in Figure 

c. Solution of the equation of motion. 
Throughout much of this section formulas of Munoz and 
Pavic[5] will be transcribed essentially unchanged, ex¬ 
cept for bringing symbols into consistency with conven¬ 
tional accelerator notation. The Munoz/Pavic formula¬ 
tion, though equivalent to various other formalisms de¬ 
scribing relativistic Coulomb orbits, is especially appro¬ 
priate for our relativistic accelerator application. Munoz 
and Pavic show that the “generalized”-Hamilton vector 


h. = hj-v + hg 9 


( 4 ) 


is especially useful for describing 2D relativistic Kepler 
orbits. Our 3D accelerator application can be formulated 
in such a way as to use only such 2D orbits. Except for 
scale factors, and a q-dependent offset, the pair {hj.,hg) 
will reduce to the conventional phase space coordinates 

(x,x')Q 


^ The overhead tildes in Eq. have been added for our later con¬ 
venience, when, for our accelerator application, we will rescale 
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FIG. 1. The bold curve shows a particle orbit passing 
through a spherical, m = 1, electrostatic bending element. 
The shaded surfaces are electrodes and the hgure is grossly 
distorted. The “Q” shown at the origin is the “effective point 
charge” that would give the same idealized electric field as 
the electrodes. 


In accelerator context the evolution of h for any par¬ 
ticular particle will be interpretable as the phase space 
evolution of that particle. We will refer to h as the “MP- 
vector” since, as far as we know, Munoz and Pavic into- 
duced it. In the non-relativistic regime h is related to the 
(nonrelativistically conserved) Laplace-Runge-Lenz vec¬ 
tor. As generalized by Munoz and Pavic, though not 
quite conserved in the relativistic regime, h satisfies a 
very simple equation of motion, whose exact solutions 
are sinusoids, even at arbitrarily large “betatron” ampli¬ 
tudes. 

In geometric mechanics jargon [7], the lattice is there¬ 
fore “integrable”. The h formalism then makes it pos¬ 
sible to define an implicit transfer matrix (where “im¬ 
plicit” implies the matrix elements depend on each par¬ 
ticle’s own conserved orbit parameters, a.k.a. “orbit el¬ 
ements”). This has the virtue of hiding the strong cou¬ 
pling between kinetic energy and horizontal position in 
an electric ring. 


the generalized-Hamilton vector slightly, for example making it 
dimensionless rather than a velocity. Dividing by a constant 
characteristic velocity will bring the formalism more into con¬ 
formity with standard accelerator terminology. Except for these 
overhead tildes, in this section, quantities and formulas will be 
largely copied from M &: P; this includes using the abbreviation 
f for df /d9. When the formulas have finally been interpreted in 
accelerator context, the usual accelerator definition of the prime 
symbol, with f standing for df/ds, where s is arc length, will 
be adopted. On the design orbit of radius rg, where 9 and s 
are both defined, one has s = ro9, and this rescaling merely 
introduces constant factors ro into the formulas. 


FIG. 2. The bold curve shows a particle orbit passing 
through a curved-planar cylindrical, m = 0, electrostatic 
bending element. The electrode spacing is g and the design 
orbit is centered between the electrodes in the y=0 plane. 


Though h is not conserved in general, Munoz and Pavic 
show (and it will soon become obvious) that h is con¬ 
served if and only if the orbit is circular. For the proton 
EDM lattice the design orbit is circular within bends. 
Also, for weak-focusing lattices, neglecting the effects 
of lattice quadrupoles situated between the bends, off- 
momentum closed orbits are also almost circlular. 

But h also has the numerical disadvantage mentioned 
earlier—transverse displacements (the normally signifi¬ 
cant aspect of accelerator dynamics) need, at least super¬ 
ficially, to be described by expressions having large can¬ 
cellations. This makes it essential to use series approxi¬ 
mations very carefully, even when the linearized paraxial 
approximation is fully valid. It is important to organize 
the evolution formulas in such a way as to exploit can¬ 
cellations analytically, for example arranging for leading 
terms to cancel analytically, rather than relying on their 
numerical cancellation. This can usually be accomplished 
(without sacrihcing exact expressions) by simplifying the 
sum of the approximately cancelling terms and keeping 
track and restoring any discrepancy. This recovers nu¬ 
merical precision for small amplitudes comparable to that 
achieved in conventional Courant-Snyder formulation. 

Another inconvenient aspect of the Muhoz/Pavic for¬ 
malism is the fact that the design, central orbit is de¬ 
generate in the sense that it is, in fact, the arc of a per¬ 
fect circle. This is precisely the condition in which the 
perihelion (or aphelion) angular orbit element is unde¬ 
fined because the radius vector length is independent of 
angular position. This ambiguity is most important in 
the time of flight calculations needed for accurate model¬ 
ing of synchrotron (longitudinal) oscillations—a calcula¬ 
tion that strains the numerical accuracy because all path 
lengths are so nearly equal. 
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These same considerations complicate the derivation of 
explicit (meaning the elements are independent of parti¬ 
cle coordinates) transfer matrices and their generaliza¬ 
tion to truncated power series transfer maps. This also 
is only a inconvenience that can be worked around. 

With qualitative discussion complete, the analytic for¬ 
mulation can begin. The Lorentz force equation in the 
TO =1 spherical case is 


dp 

dt 



(5) 


where k is the customary MKS notation for l/(47reo) ex¬ 
cept for implicitly including also the (point) charge fac¬ 
tor. On the central orbit the centripetal force equation 
is 


{pQc/e)Po 

ro 


also 

TT = ^0 


( 6 ) 


where k, like Eq, the on-axis electric field, is defined to 
be positive. In UAL, as in the MAD simulation code, 
energies like poc are always expressed numerically in GeV 
units, such that their numerical value in the computer is 
10“®poc/e, giving them GV units—an MKS unit. From 
Eq.(§, fc/e=(poc/e)/?oro and k is always expressed in the 
computer as 10“®fc/e, which gives it GV-m units—also 
an MKS unit. The on-axis electric field Eq is measured 
in GV/m units. In these units k = EqTq, numerically, 
with ro measured in meters. 

A particle’s angular momentum is 


L = r X p. 


(7) 


In terms of the design momentum po 3 'iid the design ra¬ 
dius ro, To = and 


— = Po (= 0.598379 for pEDM). ( 8 ) 

LqC 

Like k, internally (i.e. in the computer), Lqc is expressed 
numerically as 10“®Loc/e. For internal numerical val¬ 
ues to be most easily substituted into an analytic exter¬ 
nal documentation formula, factors should be grouped as 
pc/e, Lc/e, k/e, etc. or, as here, k/{Lc), and factors of 
10® (for eV to GeV conversion) put in “by hand”. Eq. ([^ 
is also useful in the form 


element, where it depends on the local electric potential) 
and the angular momentum 

L = 7^mpr^0, (11) 


are constants of the motion. It has been necessary to 
distinguish between the Munoz potential energy bVm and 
our potential energy eV, because our potential vanishes 
on the design orbit, while theirs vanishes at infinity. For 
the design orbit, using Eqs. and ([l 0 |, one obtains 


Tm,o = 


7o 


( 12 ) 


(The curious presence of 70 in the denominator follows 
from the definition of potential energy.) Using the obvi¬ 
ous unit vector relations (valid for motion with angle 0 
increasing) 


0' = —f, and f' = 0, (13) 


where primes stand for d/d9, the Lorentz equation can 
be re-expressed as 


dp k dO 
dt dO 


(14) 


Defining covariant velocity u = 7 ^v, 
to express d6/dt in terms of L, the u 
is 


and using Eq. (11) 
equation of motion 


du dO 

dO L dO 


Munoz and Pavic then introduce the generalized Hamil¬ 
ton vector 



which we refer to as the MP-vector. It is specially tailored 
so that the differential equation for h reduces to 


dh k ^ 

dO L dO 


(17) 


(The previously mentioned constancy of h on circular 
orbits is evident from this equation, since 7 ^ is constant 
on, and only on, circular orbits.) Using 


k = ^ompC^roPo- (9) 

In SXF lattice description files (discussed further in the 
accompanying paper [5]) the RF parameter V actually 
stands for 10 “® times the maximum energy change in 
eV of charge e passing through the RF cavity. 

One shows easily that both the total energy 

£ = (j^TTipC^ = £m + —, (10) 

V ry ro ro 


L ~ L kj^ also k^L k £m 

^9 5 ^9 j 2 ’ 

rupr rripr L rripr L rripC^ 

(18) 

where the new parameter k (soon to be identified as hor¬ 
izontal betatron “tune”) is defined by 



On the design orbit 


(where the “I” superscript on 7 ^ re-emphasizes that this 
relativistic factor has to be evaluated “inside” the bend 



— (= 0.801213 for pEDM). (20) 

7o 
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(Of course kq will not be controlled to this accuracy in 
practice and, in any case, quadrupoles and drifts in the 
lattice will alter the ring tunes.) Solving Eq. (18) for r, 
the orbit equation in bend elements can be expressed in 
terms of /ig; 


A 

1 + e/ig 


( 21 ) 


where 


A = 


k£ 


and e = 


M 


rrinC 


'M 


( 22 ) 


The parameter A is especially important in the storage 
ring context, because the second term in the denomina¬ 
tor of Eq. (21) is both small compared to the first, and 
oscillatory, ^s a result A can be thought of loosely as 
an average radial position. (As earlier, the “constants of 
motion” L, k, e, and A are to be referred to as “orbit 
elements”. Though constant along any single particle or¬ 
bit, they are different (though only very slightly) for each 
individual particle.) 

In component form Eq. 0 is 


matching to the known initial conditions. These new 
“orbit elements”, C and Oq, are expressible in terms of 
the previously defined orbit elements. They are now to be 
used to parameterize the true, particle-by-particle orbits 
through the bend element. (As explained elsewhere, the 
effects of the small deviations of the actual electric field 
from the Coulomb field are to be corrected for by virtual 
delta function kicks between bend slices.) Concentrating 
on hg, and defining 6 to be zero at the bend entrance, we 
have 


/io = he\e=D = C cosk6>o, 

/iq = h'g\g^Q = CksuikOq. 
These equations determine C and to satisfy 



and 


(28) 


h' 

tanK^o = (29) 

Khg 


The initial conditions for hg need to be determined from 
the known coordinates of the particle just as it enters the 
bend region. From its defining Eq. (16), 


— hg — 0, hg hj- 



(23) 



(30) 


The relativistic factor 7 ^ (equal to {£m + k/r)/[nripC^)) 
can be expressed in terms of hg. First, eliminate r in 
favor of Uq using the first of Eqs. (18), then eliminate Uq 
if favor of hg, and, finally, solve for 7 ^; 


ryl = 


£m 


k^Lc^ 


hg. 


(24) 


Note that 7 , like L, is a running, particle-specific probe 
variable, updated on every entry to or exit from a bend 
element. When inside the bend it is equal to 7 '^, but 7 
varies discontinuously as the particle passes bend edges 
(off-axis). In Eq. (18), for any particular particle, r is the 
only factor which is not a constant of the motion. We 
therefore have 


Needed for simplifying Eqs. ( 

23 

), differentiating this 

dh 


2 L dr 


2 L dx 

equation yields 



d9 

0=0 

rUpr'^ d9 

61=0 

9 J 

rupr^ as 


7 ^' = 




k2Lc2 ® 


(Heuristically, this expresses the dependence of mechan¬ 
ical energy accompanying the change in potential energy 
associated with change in particle position.) After sev¬ 
eral lines of algebra, Eqs. (23) then reduce to 


h^ — hg, 

hg = -K^hr- 


(26) 


=0 

(31) 


(25) The second of Eqs. (29) can then be used to fix 9^, but 


this requires resolving the multiple-valued nature of the 
inverse tangent function. One hopes this can be done 
by requiring r to be “smooth” (continuous, with contin¬ 
uous slope) across the bend entrance edge. Even this is 
tricky since, for non-normal incidence, there is a refrac¬ 
tive deflectior0 at the boundary. This deflection is small, 
especially near normal incidence, which is always true 


These are the equations that justify having introduced 
the generalized Hamilton vector. Their general solution 
can be written as 

hg = C cos k{9 — 9q) 

hr = ^ sin k{9 — 9q) . (27) 

where d is a running angle in the interior of the bend 
and 9o is an angle to be determined, along with C, by 


2 A “refractive deflection” is mentioned here and at other places 
in the paper. With “hard edge” bending elements the electric 
potential changes discontinuously from just outside to just in¬ 
side the edges of the element. The mechanical energy has a 
corresponding discontinuity and, as a result, so also does the 
magnitude of the momentum (which is predominantly longitu¬ 
dinal). But, since there is no transverse force, the transverse 
momentum is continuous across the boundary. Taken together, 
the result is an angular (Snell’s law like) refractive angular de¬ 
flection of each orbit at each field edge. Similar deflections occur 
at quadrupole edges. But, as well as being far weaker, the en¬ 
try and exit deflections tend to cancel in all thin elements. In 
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in our case. To resolve the inverse tangent ambiguity 
for non-normal incidence it would usually be sufficient to 
choose the value for which the slope is most nearly con¬ 
tinuous. But to be safe one has to account explicitly for 
the refractive kink. 


From Eq. (21) the radial coordinate r is given by 


A 


r = 


1 -I- ecos k{0 — 9q) ’ 


where e = Ce. 


(32) 


This is the orbit equation previously introduced in 
Eq. ([^. Borrowing terminology from solar planetary or¬ 
bits about the sun, from its structure, one sees that Oq 
is either the angle at perihelion, or aphelion, depending 
on the sign of e. In multi-million turn tracking, there are 
many opportunities to misidentify aphelion and perihe¬ 
lion or to obtain the wrong sign for e —this complication 
is associated with the degeneracy of the central orbit that 
was mentioned earlier. 


Using Eqs. (22), (29), and (32) 


= C 



(33) 


propagation plane to be the plane containing the incident 
velocity vector. Instead of keeping track of y we have to 
keep track of the orbit plane or, equivalently, the normal 
to the orbit plane. This vector changes discontinuously 
(but only by a tiny angle) in passing through thin mul¬ 
tipole elements. It is undefined in drift regions and is 
constant in bend elements. Keeping track of the orbit 
plane is covered in later sections. 

d. Rescaling of the MP-vector and updating 
the horizontal slope. For updating the horizontal 
slope component, we introduce the previously-announced 
rescaling of the MP-vector; 


h = 




(35) 


Dividing by the design velocity vq has rendered h dimen¬ 
sionless, and the further factor of 70 in the denominator 
simplifies matching to conventional lattice function rep¬ 
resentation. At this time we also switch definitions of 
the prime operator notation so that f' = df/ds, where 
s is arc length along the design orbit. With this revised 
notation Eqs. ([26| become 


By expressing the small parameter e, proportional to the 
small Hamilton vector value and slope (or rather their 
quadrature sum), this formula is essential to the imple¬ 
mentation of the Muhoz/Pavic approach. (These “small” 
quantities can even be infinitesimal in the analytic treat¬ 
ment.) 

From these equations one has obtained simple analytic 
parameterizations for x{9) and 7^(0). In particle tracking 
one knows to high accuracy and wishes to find Xoutj 
also to high accuracy. Their difference is given by 


^out ^ii 


A A 

-Ir-^0-IW + ^0 

1 + e 1 + e /iq 

ho — hg 

-- A€. 

(1 -|- e hg){l + e ho) 


(34) 


This manipulation has suppressed the harmful cancel¬ 
lation exhibited in Eq. (32) which gives r rather that 
the radial displacement from the design orbit. Eventu¬ 
ally, when deviation from I/r^ electric field dependence 
is modeled by thin virtual quadrupoles, to improve pre¬ 
cision, one will slice the bends finer and finer. 

In the context of accelerator physics one can ask for 
the similar evolution equations for vertical coordinate y. 
But this question is inappropriate; the Kepler orbit lies, 
by definition, in a single plane and we always choose the 


ETEAPOT the refractive corrections are included at all bend 
edges and neglected at all other edges. At bend field edges there 
are also fringe field corrections (which are especially important 
sources of spin decoherence.) The fringe field and refractive cor¬ 
rections are evaluated separately and simply superimposed at 
every bend edge. 


hL 


h'. 


— hg, 

ro 



(36) 


Einally it is possible to make contact with more familiar 
accelerator variables, such as Px and Ps, the radial and 
longitudinal particle momentum compon ents . 

From the definition of hr in Eq. (16), after re¬ 


scaling (35), the transverse momentum is given by 


dx 


Px 


Px = lompVo hr, or x' = — = p[l] = — = hr, (37) 


Vz 


where po = 'joixipVo is the design momentum and p[l] 
is the symbol used for x' internally in the UAL code. 
(Some components in the chain of equalities in Eq. (37) 
may assume paraxial approximation, even if for no reason 
other than that the definition of x is unambiguous only 
in the paraxial limit.) So hr is nothing other than the 
phase space coordinate conjugate to x, namely dx/ds. 

When applied at the exit of a bend, the value of x' 
given by Eq. (1^ would better have been called _, 
since the refractive compensation associated with exiting 
the bend would still need to be made to produce p[l] 
valid just outside the bend exit. 

Potential loss of numerical precision is always an issue. 
A compact, numerically accurate way to update the 2D 
phase space coordinates is to work directly from Eq. (32). 


Eor bend angle A9, the argument 9 — 9o is A9. The 
increment in x from input to output is 


^out ^in — Ac 


COS k9q — COS kA9 


(1 -I- ecos kA9){1 + ecos k9o) 


(38) 
























Differentiating Eq. (321 produces 


— {9) - \c:z smK{e-eo) 

(l + £ cos k{9 — 9o)Y 


(39) 


Including the change of independent variable from 9 to 
z, the increment from input to output is 


dx 

dz 


out 

Ae/c 


dx 

dz 


sin kA 0 


sin k9o 


vq \{1 + ecos kA9)^ (1 + ecos k 0 o))' 


(40) 


These formulas avoid taking the difference of nearly equal 
terms. 

e. Pseudoharmonic description of the motion. 

In general the MP plane will be close to, but not exactly 
identical with the horizontal lattice design plane. In this 
section, for simplicity in making contact with Twiss func¬ 
tions, these two planes will be treated as equivalent. 

At this point we wish to correlate the dynamic quanti¬ 
ties introduced so far with more familiar (to accelerator 
scientists) quantities such as Twiss functions, betatron 
phase advances, transfer maps, and so on. The rela¬ 
tionships are especially simple for a combined function, 
weak-focusing ring; this will be described before develop¬ 
ing the full Twiss function formulation needed for sepa¬ 
rated function lattices. In general, the definition of lat¬ 
tice functions within electric bend elements will be more 
complicated, but also not very important for the EDM 
experiment. 


Eq. (36) yields 


In practice this issue is largely academic. In any real¬ 
istic ring there are many drift spaces, more or less uni¬ 
formly distributed. Especially with the weak focusing ex¬ 
pected in EDM rings, the ranges of /3 function variations 
will be quite limited. Simply interpolating the (slightly 
ambiguous) beta functions in the bend regions from their 
(reliably known) values in the drifts should be sufficient 
for most purposes. 

As already stated, our electric lattice model ad¬ 
mits only uniform, inverse square law bending elements 
(though with artificial thin trim elements). Within such 
a bend we approximate /?' = d/3/ds = 0 and horizon¬ 
tal “betatron” oscillations of amplitude a are described 
by 

X = costp, 

x' = —sin sin -ip. (43) 


Differentiating once more yields 

x' = —costp = —P~‘^x. 


(44) 


(45) 


Substituting this into Eq. (36) produces 

X = - he. 

ro 

Just as hr and he vary “in quadrature”, so also do x and 
x'. The ratios of their maxima are 

hf) rr\pt.yc , ^max ^9 ,] 


^0 1 ^max 

h' ~ ~k’ ~x’ ~ 

"-e.max ^ -^max 

Collecting results, we have 

P = —, X = -X he, 


^ 0 ,max 


(46) 


and X — hr. 


(47) 


he 



(41) 


In traditional Courant-Snyder (CS) formalism the beta¬ 
tron phase Ip and the positive-definite lattice function 
/3(s) are related by 


^ 1 

ds /3(s)’ 


(42) 


Even with 9 = s/tq, because /3(s) is not necessarily con¬ 
stant, the angles ip and 9 = s/rg, though monotonically 
related, are not strictly proportional. The virtue of the 
beta function formalism is primarily due to its ability to 
describe each orbits in a complicated lattice as a sinu¬ 
soidal function of unambiguous phase advance ip. 

Shortly we will have to admit that our generalization 
of the CS formalism to electric lattices will be limited, in 
principle, to describing orbits outside bending elements, 
where the electric potential vanishes. The parameter /3 
introduced now will be applied inside bends, and it is 
assumed to be independent of s. Furthermore it will not 
match up smoothly with the adjacent outside P func¬ 
tions. This reflects the small discontinuity in particle 
kinetic energy at bend edges. 


(Remember that these x and x' values are only the pre¬ 
cise Frenet-Serret laboratory coordinates when the M-P 
plane coincides with the horizontal laboratory reference 
frame, which is always a good approximation.) To sum¬ 
marize, one sees that, to linearized approximation, the 
components of the MP-vector are, except for scale fac¬ 
tors, identical to the betatron phase space components. 

The parameter P has been used in this section only 
to produce familiar-looking formulas. For a weak focus¬ 
ing ring we now see that this P need never have been 
introduced, since it is just an abbreviation for tq/k. It 
does not depend explicitly on s, but through its k factor 
it depends on 7 , x and x', and is therefore different for 
different particles. 

/. Revolution period. For modeling longitudinal 
dynamics it is critical to obtain the revolution period 


Trev. to good accuracy for every particle. From Eqs. (11 1 , 


(24), (27), and (32) one has 


dt 

dB 


£mX^ 


Lk c _|_ g( 3 Qg ^( 5 ) _ 

mpkOX^ cos k{9 — 9q) 

( 1 -|- e cos k{ 9 — 9o)j 


( 48 ) 
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To obtain time of flight dt/dO has to be integrated over 0. 
Both integrals can be evaluated in closed form to obtain 
t as a function of 9. However the formulas are compli¬ 
cated and their evaluation is numerically treacherous, for 
reasons discussed below. 

For the study of longitudinal dynamics and syn¬ 
chrotron oscillations what is needed is the deviation of 
the time of flight from the time of flight of the design par¬ 
ticle. This is a very small difference of two large numbers. 
Eq. (48) can be recast as 


t — to — 4l 


dA9 

{1 + e cos kA9)'^ ° 

cos kA 9 dAO 


BC 


dAO- 


(1 -I- ecosKAd)2' 


(49) 


Here Ad = d — dp, the factors A and B (both of 
which deviate little from constancy and can never be¬ 
come “small”) replace the other factors in Eq. ( [4^ , and 
Aq is the values of A on the design orbit. All of the fac¬ 
tors A, H, Aq, e, and k are constant on any particular 
orbit being followed, and all have values very close to 
their values on the design orbit. The only variable factor 
is Ad. The first two terms cancel approximately, which 
makes the evaluation of the first term critical. The third 
term requires no special treatment since the C factor is 
already differentially small. The cancellation tendency 
can be ameliorated by adding and subtracting a term 
Aq / dA6 = AqA9 to produce 


t-to=e 


C (—2A-I-H/e) cos/cAd — Aecos^ kAO 

I (1-I-ecosKAd)^ 

+ (A-Ao)Ad. 


dA0+ 

(50) 


Here C = e/e has been obtained using Eq. (32). In 


Eq. (50) all terms are differentially small and the small 


differences of large numbers have been moved outside 
the integrals. For the proton EDM experiment, since 
|e| never exceeds 1/1000, the integrand can be expanded 
in powers of e before integrating. A typical leading term 
in these integrals is Ad multiplied by a constant factor. 
While using this time of flight formalism ETEAPOT in¬ 
cludes this term as well as terms with extra powers of 
e up to and does not check automatically that this 
last term is, itself, negligible. An alternate, seemingly 
more robust, hybrid time of flight calculation, now used 
by default in ETEAPOT, is desribed later. 

A subtle feature of the time of flight evaluation results 
from the fact that the ellipticity e vanishes on perfectly 
circular orbits (of which the design orbit is one). This 
causes the angle do, which is the angle to perigee, to be 
indeterminate because perigee is indeterminate on a cir¬ 
cular orbit. Eq. (32), the second of Eqs. (22) and Eq. (29) 


force e to be non-negative. And yet the major and minor 
axes can switch from prolate to oblate when an arbi¬ 
trarily small kink is applied to an orbit for which e is 
close enough to zero. When this happens the angle do 
advances discontinuously through tt/2. When this hap¬ 


pens the angle d also advances discontinuously through 
the same angle. 

These changes do not affect the integrals in Eq. (50) 
which depend only on d—dg. They do, however, make the 
time of flight calculation hyper-sensitive to the determi¬ 
nation of do, which can change erratically in progressing 
from one bend element to the next, after passing through 
a straight section which possibly contains a quadrupole. 
A discontinuous change in dg can change the analytic, 
closed-form indefinite integrals used in applying Eq. (481. 


IV. TRANSFORMATION FROM LOCAL 
FRAME TO MP FRAME. 

g. Determination of the MP bend plane. In 

order to reserve x for radial displacement in the Muhoz- 
Pavic (MP) bend plane we now use overhead bars to des¬ 
ignate conventional Courant-Snyder (CS) coordinates. 
At a location where the longitudinal coordinate has been 
shifted to be z = 0, the particle position {x, y, 0) is given 
by 



Scaled momentum components pj = {px/po)/i^ + Sp/po), 
and Py = (j)y/po)/{l -I- Sp/po), are defined by 

"2 

p =P 2 X-hpj,y-f (l - y - z, (52) 

where terms of fourth order and beyond have been 
dropped. (For the EDM experiment such terms have 
numerical values such as (0.01/20)^ « 10“^^.) This has 
been arranged so that p is a unit vector. The p[l] and 
p[3] MAD/UAL coordinates being tracked in the com¬ 
puter are close to, but not exactly equal to px and py. In 
fact 


Px = 


P[l] 


I -I- Sp/pQ ’ I -I- 5p/po' 

Scaled angular momentum components are then given by 


P[3] 


(53) 


L =-=-^ = det I I -I- x/tq y/r^ 0 

roPc ropQ 


Px 


]L 

ro 


Pzi- (l + ^'^Pzy + (py + 


Py Pi/ 

X ^ y ~ \ ^ 

— Py -Ps Z- 

ro ro 


(54) 


The scaled central angular momentum is Lq = —y, and 
the perpendicular component Lj^ is 


L, = 


ro 


'l_§L_Ph^ 


)-( 


X ^ y 

Py-\ - Py -P$ I z 

ro ro 


= Sin (z >2 X — sm 0sz, 


(55) 
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which defines two (small) angles the tilted bend plane 
makes with true vertical. This equation determines the 
(small) rotation angles about transverse axes that rotate 
the horizontal plane into the MP plane. The trigonom¬ 
etry is not exact. But it is accurate enough for the ex¬ 
tremely small angles involved. 

h. Merging 2D MP-tracking into 3D. With the 
orbit geometry being worked out in the 2D MP plane, it 
remains necessary to transform back and forth to the 
local 3D Frenet frame. The goal of a tracking calculation 
is to obtain the output particle radius vector Tout for a 
given input particle radius vector Tin. The corresponding 
unit vector evolution is from fin = {x-m, yin, 0) to font- We 
calculate this unit vector evolution first. 

On the design orbit the input and output vectors are 
din = (1,0,0) and dout = (cos 0,0, sin 0), both lying in 
the horizontal (x 2 ;)-plane. With the design bend angle 
defined to be 0, one has 

din • dout = cos 0. (56) 

The particle-specific bend angle (j), in its own bend plane, 
satisfies 

f in • fout = cos (/), (57) 


V. TRANSFORMED DEPENDENT VARIABLE i 
ORBIT REPRESENTATION 

i. Electric sector bends with field index m. 

The Munoz and Pavic, Hamilton vector approach de¬ 
scribed so far, though ideal for fast exact particle track¬ 
ing, is less well matched to developing the Courant- 
Snyder Twiss function description conventional in accel¬ 
erator theory. A more conventional approach is taken 
in this section. As well as enabling a truncated power 
series formalism for electric rings, the results can be 
used for corroborating results already obtained. This in¬ 
cludes independently confirming the integrable dynam¬ 
ics demonstrated for to= 1. Though discovered using 
the Munoz/Pavic approach, this remarkable integrabil- 
ity property must necessarily be derivable also in the 
paraxial-based approach. Initially this treatment will be 
generalized to include field indices other than m = 1. 

To produce weak vertical focusing the radial electric 
field has to fall off as with m > 0. The m = 0 

case is singular, and leads to a logarithmic potentialj^ 
Introducing design radius tq and central field —Eq the 
electric field for y=Q is 


where (j) is close to, but not equal to, 0 in general. Letting 
a and b = ±|6| = be the components of font 

along dout, one has 


fout = (a COS0, i-s/l — a^, asin0). 


(58) 


(For clockwise orbit) the unit vector normal to the bend 
plane. 


n = —- 


(59) 


is known from calculations based on the particle coordi¬ 
nates at the entrance to the bend element. To determine 
a we require 

0 = fiin • font = rixa cos 9 Euy^/l — + n^a sin0. (60) 

Solving for a, 

Uy 

a = 

\ cos^ 0 -f 2nxnz cos 0 sin 0 -|- n1 sin^ 0 

(61) 

The angular advance (j) then satisfies 

cos 0 = fin • fout = Xina COS 0 -|- ijinb. (62) 


The radial coordinate is then given by Eq. (21) which is 
repeated here; 

^ (63) 


1 -|- chg 


where hg is given by Eq. (27). The output radius vector 
is then given by 

Tout = J'fout- (64) 


In this way the evolution can be completed using just 
Courant Snyder coordinates. 


E(DO) = -Eo^f. (66) 

The electric potential P(r), adjusted to vanish at r = rg 
was given earlier in Eq. (§. The independent (longi¬ 
tudinal) coordinate s is to be replaced by the angular 
coordinate 0 

0=—. (67) 

^•o 

Following standard treatments of relativistic Kepler 
orbits[8][9][IO]|II], we change dependent variable from 
a:(s) = r — tq (with independent variable s) to a (di¬ 
mensionless) dependent variable ^(0) (with independent 
variable 0); 


^ a; X 

r Tf) + X 
To dr 






( 68 ) 


Note that ^ is proportional to x for small x and that, 
notationally, x' = dx/ds and = d^/dO. The present 


^ Integrating the electric field from r = lto‘r = l + A one recon¬ 
structs the potential from the electric field—call it V. Its value 
at r = 1 -h A is 

(1 + A)-^ - 1 

V{1 + A) = ^ ^ ^ - = ln(l - A) -h 0(m). (65) 

m 

The logarithmic potential can be seen to be a degenerate form 
required for m = 0, where E oc 1/r. The electric field corre¬ 
sponding to this potential is oc This shows how V 

approaches the logarithmic potential in the limit of small m. 
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discussion is limited to planar orbits, in which case the 
definition of cc by r = rg + x is the usual Frenet-Serret 
definition. For 3D motion this definition is not quite ex¬ 
act but realistic vertical amplitudes will always be small 
enough that the effect on x of projection onto the hori¬ 
zontal plane will be neglected. The identity 


For TO=0 the invariance of the system to translation along 
the y-axis, causes py to be conserved, and to invariance 
under rotation around the central axis, causes Ly, the 
vertical component of angular momentum, to be con¬ 
served. For a particle in the horizontal plane containing 
the origin the angular momentum vector is 


— = 1 - ^ (69) 

r 

will be prominent in subsequent formulas. Inverse rela¬ 
tions are 


rxp = 


rupry 


:d+- 


The design orbit angular momentum is 


y = Le6+Lyy. 

(75) 


X = 


^ 0 ^ 

1 -r 


and x' 


rod^/ds 


(i-c)^' 


(70) 


With ^ regarded as a function of 9 and a; as a function 
of s, the abbreviations p = d^/d9 and x' = dxjds have 
been used. 

For simplicity in describing the approach to be taken, 
we temporarily specialize to the to=0 “cylindrical” case. 
The electric potential V{r), adjusted to vanish on the 
design orbit, is 

T 

V{r) = — = Sorglnr - Fiorglnro. (71) 

ro 


Lq = mpC/3oJoro- (76) 

We seek the orbit differential equation giving dependent 
variable r as a function of independent variable 9. The 
(conserved) total proton energy £ is the sum of the me¬ 
chanical energy and the potential energy 

£ = ^JprP +^ 0 ^ +PyP + -\- eV(r). (77) 

Squaring this equation yields 

{£ - eV(r))^ = pIP +pIP + m-pp. (78) 


The design parameters are related by 

eEgro = (3opoc = mpp(yo - —) , (72) 

70 7 

where p, v, and /3 are proton momentum, velocity, and 
v/c. The momentum vector components are defined by 


Tp=PrV+ pe9 + py-y (73) 

_ TTipr ^ ^ mpr9 ~ rUpy , 

y/X — jp ^JX — jf? ^JX — jp ^ 

The electric force alters only the radial momentum com¬ 
ponent 


dpr 

dt 




The terms on the right hand side of this equation can be 
expressed in terms of r, dr/d9, and conserved quantities. 
From Eqs. (Fz^ and ([7^, 


P9 = —■ 
r 


(79) 


(80) 


Pr can be expressed similarly using Eq. (731; 

Pr r X dr Ly dr 

pe r9 r d9 r^ d9 

Making these substitutions yields 
(£-eE(r))'= +^+py + my. (81) 


When the field index m was first introduced it was visualized as being a small number, just large enough (and 
positive) to produce some vertical focusing. But the subsequent analysis has not relied on m being small. Generalizing 
the previous formulas for arbitrary m values, and expressing the in-plane momentum in terms of angular momentum 
Ly, ^ and p, the mechanical (total minus potential) energy can be expressed in terms of particle mass and momentum 
components. This generalizes Eq. (81) for arbitrary m; 



eEpro 

m 


+ 


eEpro 

m 


(1-er 


i^fd^\ 

[d9) 


+ 


Ly 


(1 


-y 


I 2 2 

+ PyC 


2 4 
rripC . 


2 


(82) 
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Differentiating this equation with respect to 9, and simplifying, 

^ / eEpr^ _ V\ 1 _ e'^E^r^ 1/m 

dO"^ \ mj (1 —(1 — 


1-?-1 

£ Sje Eoto I 

f Eoro \ 

'^-1 
m ) 

1 1 1 

t Eoro \ 

2 1/to 

\Lycl{ero) Lyc/{ero) 

\Lyc/{ero)) 

' (I - 0^""* 

KLyc/(ero)J 



^ L| _ Ll 1 _ L| /3g 1 

£o LI LI m) (1 - 0^-™ LI m (1 - ' 


(83) 


This equation can be reduced to quadratures but, the subsequent integral cannot be performed analytically for non¬ 
integer m values, and only with difficulty, if at all, for m values other than 1. 


j. Specialization to to = 1 Planar Orbits. The 

formulas just obtained simplify spectacularly for to=1 . 
For example, the electric potential (defined to vanish for 
^ = 0) is given by 


V{0 = Eoro^- 


(84) 


To obtain the orbit equation we can start from Eq. (831, 
specialized to the spherical electrode geometry shown in 

Fig.[T] 

Eq. ( [8^ becomes especially simple for motion in the 
y = 0 plane with to=1 ; 


d6i2 LI So 



^ = (85) 


The last form implicitly defines abbreviations Qx and ^co¬ 
lt also introduces a representation of transverse displace¬ 
ments as the superposition of a “betatron” part and an 
“off-energy” part. 

It can hardly be over-emphasized how simple this equa¬ 
tion of motion is. When expressed in terms of ^ the mo¬ 
tion is not only simple harmonic, it is simple harmonic for 
arbitrarily large oscillation amplitudes. The only reser¬ 
vation is that ^ cannot exceed 1, at which point r=oo. 

For our application the coefficient is positive, which 
means there is horizontal focusing with horizontal tune 
Qx given by 


Qx 



( 86 ) 


The off-energy central orbit is given by 


^co 


1 



LI £o)^ 


C = 0 . 


(87) 


With Ly and S allowed to vary, (especially in RF cavities, 
but not within electric bend elements) the parameters Qx 
and Xco have to be regarded as locally, rather than glob¬ 
ally, defined. It is important to notice though, that the 
parameters entering the definitions of both Qx and ^co 
are invariants of any given particle’s motion. The quan¬ 
tities Eq and /do are obviously invariant because they are 


properies of the on-energy central orbit. The quantities 
£ and are invariant by conservation of energy and 
angular momentum, but they are, in general, (slightly) 
different for different particles. 

k. Particle-specific transfer matrices. In prac¬ 
tice there will be drift regions, quadrupoles, RF cavities 
and other apparatus in the ring, not to mention artihcial 
thin quadrupoles inserted within the bend to compensate 
for actual field index differing from to= 1. Eq. ( [85| will 
be applied only within individual sector bend elements. 
But the simplest possible application of the formulas has 
a single element bending through 27r—which is to say, 
making up the whole ring. The x motion described by 
linearized treatment is only sinusoidal for small ampli¬ 
tudes. But, since the denominator is a periodic function 
of 0, x{9) has the same period. For the frozen spin value 
jdx = 0.59838, with orbits close enough to circular that 
L « Lqj the tune is Qx = 0.8012. (Notice that the 
parameter Q and the Munoz parameter k are identical 
parameters.) 

Our equations, being fully relativistic, are valid in 
the nonrelativistic limit ,5o —t 0. In the nonrelativistic 
limit the orbits are simply Keplerian planetary orbits, 
which are known to form closed ellipses, corresponding 
to Qx=l- So one can sot that the shift of Qx away from 
1 is a relativistic effect o 

An orbit having initial vertical displacement or slope 
leaves the horizontal (design) plane. It nevertheless lies 


^ As a curiosity, one can evaluate the corresponding tune shift of 
the orbit of planet Mercury around the sun. We pretend the or¬ 
bit is circular (when in fact its eccentricity is e = 0.21, yielding, 
for fixed total energy, (L/Lo)^ ^ 1/(1 — ri) fa 1.04.) Mercury’s 
mean orbital velocity is 47.87km/s so do = 1-596 X 10“^ and 
Qx ~ 1 — do/2 = 1.27 X 10“®; this represents a negative frac¬ 
tional advance of the perihelion each turn. (Our electrostatic 
formulation amounts to accounting for the relativistic increase 
in inertial mass with no corresponding increase in gravitational 
attraction.) Since its orbital period of 0.241 years. Mercury 
completes 100/0.241=415 revolutions in 100 years. Meanwhile, 
its Einsteinian general relativistic perihelion advance is 43 arc- 
sec. Expressed as a fractional deviation per revolution, this is 
43/(360 X 3600 x 415) = 7.99 x 10“®. It appears, therefore, 
that the perihelion advance is reduced by about 16% by the spe¬ 
cial relativistic inertial increase. This is roughly consistent with 
Chandler ll2l . 
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in a single plane, and its evolution is the same as if it were 
in the design plane. From this point of view Qy = 0. On 
the other hand, if one insists on interpreting non-zero y 
values as vertical betatron oscillations, the vertical tune 
is Qy = 1. This is an unusual example of “tune aliasing”. 

For pure horizontal betatron oscillations 

Q = Q, = (88) 

For propagating the radial displacement ~ Cco, 

one can introduce cosine-like trajectory C^{9) satisfying 
C^{0) = 1, C^(0) = 0 and sine-like trajectory S^{9) sat¬ 
isfying 54 ( 0 ) = 0, 5^(0) = 1. They are given by 


C^{9) = cos{Q9) 

C'^{9) = -Q sin{Q9) 

S'^{9) =cos{Q9). (89) 


For describing evolution of (^, ^') from its initial values 
(Ciii) ■Cin) at 0 = 0 to its values at 9 one can use the “trans¬ 
fer matrix” defined by 


Me (0) 


/Q(0) S^i9)\ 

1,^(0) S'^{9)) 


( C S/Q\ 
\-SQ C ) ’ 



where C = cos(Q9) and S = sm(Q9), to give 

or, spelled out in inhomogeneous form, 

m = Cco +c (6„ - Cco) + (s/Q) c(„, 

^'(0) = -5Q(Ci„ - Cco) + C^'n- (92) 


The subscript ^ on Mj and its components is to serve 
as a reminder that this matrix is particle-specific, and is 
therefore not a transfer matrix in the conventional, lin¬ 
earized, particle-independent, sense. In spite of appear¬ 
ing superficially to be linearized, describes nonlinear 
motion, and is exact for all amplitudes. Eq. (921 can be 
written more compactly as 





_ / C S/Q] 
- [-SQ C ) 



(93) 


1. Propagation Through a Sector Bend Starting from known coordinates just preceding the 

entrance to a sector bend, where the electric potential vanishes, one presumably knows £ and, therefore. 


pf^_=£^/c^- 


2 2 
m^c , 


V^- 

r z,\n— 


yin— 


5 Px,in — ^ in— \/ in— 


(94) 


(In our hard edge approximation) just past the entrance, x, £, and Px are unchanged, but the other dynamic variables 
have changed. In particular, using conservation of Px, 


Pin+ = 


^£ 
V C 


Epro Xin/C N^ 
ro + Xin ) 


2 2 

- , 


Pz;m+ — Pin-\- P 


2 

a;,in’ 


^in+ 


Px,in 



(95) 


This has determined The last expression has been arranged to give not just the magnitude, but also the sign 

of x[^,. (The argument of the square root can never change sign.) Notice that this calculation has included the (very 


small, because of near-normal incidence) “refractive” bending at the interface. Using Eq. (68), initial coordinates 
(^inj ?in)i just inside the sector bend can be obtained; 


6 n = 


(Tin- 


ro -I- Xin- 


4 = 


(ro + Xin-)2 ' 


(96) 


Note also that Ly, which is proportional to (rp + x)pz changes according to 


^y,in+ 


Pi.m+ 

2 • 
PzAn- 


(97) 


with both numerator and denominator known from earlier formulas. From Eqs. (86) and (87) it can be seen that Ly 
is the only variable parameter in the definitions of Q and Xco) (not counting changes in £ occurring in RF cavities). 
Q and x^o can be updated accordingly. With the updated version of Xco one can determine ^in — x^o which is needed 
as an input for the bend element transfer matrix. 
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From the known bend angle Ad, cosine and sine factors C and S can be calculated, and all elements of transfer 
matrix determined. Then, using Eq. (93), one obtains the components (^out,^out) just before the exit from the 
sector bend. 

To get out of the bend element the steps taken at the input have to be reversed. 

I 


VI. TIME OF FLIGHT AND LONGITUDINAL 
PHASE SPACE DYNAMICS 


The treatment of orbits has been based on positions, 
slopes and momenta, with no need for velocities. To cal¬ 
culate time of flight, velocity is required. These calcula¬ 
tions will proceed in parallel with the orbit calculations 
but they are described separately here to emphasize the 
essentially different evolution of position r(d) and ct{9), 
which is the arrival time (expressed as a length) relative 
to that of the central particle. With all quadrupoles and 
multipoles treated as thin in ETEAPOT, the only con¬ 
tributions to time of flight come from drifts and electric 
bend elements. The only essential complication comes 
from the dependence of particle speed on position. 

Even apart from their different velocities, with all or¬ 
bits very nearly parallel, the path length is very nearly 
the same for all particles. Yet synchrotron oscillation 
dynamics depend sensitively of their time of flight differ¬ 
ences. Furthermore, the main reason the spin coherence 
time SCT can be long is the strong tendency for excess 
spin precession to cancel over complete synchrotron os¬ 
cillation cycles. To the extent this cancellation depends 
nonlinearly on synchrotron amplitude the SCT is im¬ 
paired. These considerations make time of flight calcula¬ 
tion difficult, yet important. 

As shown previously, in the m = 1 Muhoz/Pavic treat¬ 
ment the time of flight integrals can be evaluated ex¬ 
actly and in closed form. But the resulting formulas 
are quite complicated, and numerically treacherous, be¬ 
cause they contain (multiple-valued) inverse tangents and 
cotangents. Formulas are given in the ETEAPOT user 
manual [T5]. 

In the paraxial approach (in spite of previously having 
been shown to be “integrable” for m = 1) the time of 
flight integral cannot be expressed in closed form. But, 
after power series expansion, the integrals are elemen¬ 
tary. The resulting power series is so rapidly convergent 
that accuracy to machine precision is quick. We have 
found this approach more robust than the Munoz Pavic 
approach for analysing longitudinal motion. This formu¬ 
lation is described next. 

m. Kinematic variables within electric bend el¬ 
ements. For the m = 1, inverse square law case being 
emphasized, the electric field and electric potential are 
given by Eqs. ([^ and (|^; 

E(C) = -Ao(l-0"f, 

^(C) = Eorot (98) 


The latter equation, along with Eq. (92), permits the 
potential energy to be expressed in terms of 9 and initial 


conditions. 

V{9) = Eotq (^^co + (Cin-Cco) cosQ9+^ sinQ9.j (99) 
Then 7 ( 0 ) is obtained from 

7(0) = —^- ^°2°/ (^co+(6n-'fco) cosQ0-k% sinQ6» 

rUpC^ TUpC^je V Q 

( 100 ) 

From this formula one can obtain (3(9) using 


(3(9) = y/T^W- 


( 101 ) 


From the y-component of Eq. (75) we also have, for mo¬ 
tion in the horizontal plane. 


d9 

dt 


-L 


V 


rupr^j 


( 102 ) 


Warning: Note that, with right-handed Frenet coordi¬ 
nates, and clockwise orbit (which we now assume) Ly is 


negative. This accounts for the negative sign in Eq. (1021. 


The right hand side of this equation can now be expressed 
in terms of 9, invariants and initial conditions. The neg¬ 
ative sign {—Ly) is specific to clockwise orbits, with 9 
increasing along the orbit. There seems to be no univer¬ 
sally accepted sign convention distinguishing clockwise 
and counter-clockwise beam evolution in storage rings. 
This sign ambiguity causes the tracked time of flight vari¬ 
able ct to also be ambiguous. The appropriate RF phase 
is then, similarly, ambiguous. 

n. ^ evolution and convergence estimates. Us¬ 
ing Eq. (68), initial conditions (■Co, Co) have been ex¬ 


pressed in terms of initial x conditions; 


c = 


ro 


C' = 


raX 


(ro 


\2 ■ 


(103) 


The evolution of x is obtained using Eq. (921, whose up¬ 
per equation yields 

C = Cco + (Cin - Cco) cos sin Q9. (104) 

This formula can be used to give C as a function of 9, for 
example when calculating times of flight. Then 


x{9) = 


roi{0) 

1-C(0)’ 


and r{9) = rg -I- x{9). (105) 


For the proton EDM experiment the value of radius rg 
will be, say, 40 m. An initial value ccg = 1 m would be un¬ 
realistically large. It is better therefore to use centimeter 
units. Then rg = 4000 cm, and the initial displacement of 
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a cosine-like trajectory is 1 cm. With the electrode spac¬ 
ing being 3 cm, a unit-amplitude, cosine-like trajectory 
happens to have something like a maximal amplitude. 
Dropping terms of one higher order, for example from 
quadrupole order to sextupole order, therefore amounts 
to making an error of about one part in 4000. An er¬ 
ror of this magnitude is likely to be unimportant unless 
some resonance causes its repeated effect to accumulate 
constructively over many turns. 

o. Time of flight through bends. A formula for 


the time of flight has already been given in Eq. (48 1 . 


Since this form is integrable, it gives the exact flight time 
in closed analytic form. However it has the disadvantage 
that the angle Oq (which is the angle to “perihelion”) is 
hard to determine unambiguously, because it is multiple¬ 
valued, and because perihelion is indeterminate for cir¬ 
cular orbits. The most important circular orbit is the 
design central orbit and, depending on elements encoun¬ 
tered in the ring, the local angle to perihelion can vary 
erratically. 

It is useful, therefore, to have an alternate, highly ac¬ 
curate, though not necessarily analytically exact, form, 
to accompany the Munoz form developed earlier. From 
Eq. (102 1 the time of flight dt, expressed as a distance 


d{ct) is given by 


1 


de 

where 


LyC 


7 = 


-Lyc/e (1-C)‘ 


{£/e - AoroC). 

(106) 


? = Cco + (6n - Cco) COS Q0+^ sin QO. 


(107) 


Note that the variable factor has been replaced using 
Eq. (69) and that the parameters Q and ^coj though in¬ 


variant, are particle-dependent as in Eqs. ( 86 ) and (87). 
(The proliferation of “/e” factors is for convenience in 
converting electron-volts to MKS energy units.) For¬ 
mula (|106|) can be checked on the design orbit where 


^ = 0, using Lqc = -mpC^^oPoro] 


dcto 


1 


hompC^le) = 


ro 


/3o’ 


(108) 


0 ' 

'yo/doroiripc'^le 

which is correct. The time of flight through an element 
with angle A9 is given by 

r^S/e-EoroO 


ct = 


1 


-Lyc/e (1-0' 


dO, 


(109) 


Jo 

(where, for clockwise motion, —Ly is positive). Over the 
corresponding sector the time of flight of the design par¬ 
ticle is 

/■A 6 » 

cto = / -^de. ( 110 ) 

Jo Po 

The increment in time relative to the design particle is 

"r g(£/e-AoroO _ 1 _ 

-Lyc/e (1 - if 

( 111 ) 


ct — ctf) = 


/o 




The integrand can be expanded in a series in 0 
rl{£/e - Eor^i) 1 


ro 


-Lycje (1 - if Po 

where 
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( 112 ) 
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(113) 


All the required integrals (over the range from 0 to A0, 
which is the bend angle) are elementary—the integrands 
are integer powers of sin Q6 or cos Q9^ and the series in 
Eq. (112) is rapidly convergent. It is important to re¬ 
member that quantities, in this case Ly, which are dis¬ 
continuous in the transition from just outside to just in¬ 
side a bend element, have to be updated to the value 
inside before evaluating these coefficients. 

Because C is necessarily positive while ^ can have ei¬ 
ther sign, the “quadratic” term proportional to C may 
compete with the preceding “linear” term proportional 
to i, which can have either sign and will tend to cancel 
on the average. As has been stated repeatedly, since the 
coefficients are particle-dependent, this formula has to be 
performed individually for every particle at every bend, 
p. Time of flight through straight sections. 
Different particles also have different flight times through 
drift sections or multipole elements of length £. The path 
length excess is approximately 


e{t to)path — 






2£Po 

and the off-speed correction is 


2/3o 


C(t to)vel. 

The total effect is 

— ^o)straight = 7 ;^ ( -1 + ^ 


(114) 


(115) 


(116) 


or, more directly, 

ft - to).traight = (117) 

q. Time of flight due to vertical oscillation. 

For bends only the projection of the orbit onto the hori¬ 
zontal plane has been found so far. It remains to evaluate 
the time of flight ascribable to vertical betatron motion. 
We assume there is no direct coupling between x and y. 
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We also neglect the effect of the small speed variations ac¬ 
companying horizontal oscillations. Also neglecting the 
small bending in sufhciently-finely sliced bend elements, 
we treat them as straight. With these approximations 
the bend element is treated as a drift as far as vertical 
oscillations are concerned. The off-velocity effect has al¬ 
ready been included in the horizontal calculation. Copy¬ 
ing from Eq. (116), 


tQ)vert. 


2/3o ■ 


(118) 


VII. LUMPED CORRECTION FOR FIELD 
INDEX DEVIATION 


Analytic propagation formulas have been given for 
1/r^ radial dependence of the electric field. The ac¬ 
tual radial dependence will deviate from this, with to, 
the deviant held index, having been dehned in Eq. ([^. 
E ^ l/r^’*'™, so TO = 1 for the Coulomb’s law field depen¬ 
dence. For this, the “spherical” case, the orbit equations 
have been solved in closed form. 

The ETEAPOT strategy, even for to 7 ^: 1, is to treat 
sector bends as thick elements with orbits given by the 
analytic to = 1 formulas. The m ^ 1 case is handled by 
inserting zero thickness “artihcial quadrupoles” of ap¬ 
propriate strength. As with TEAPOT, this approxima¬ 
tion becomes arbitrarily good in the thin slicing limit. 
Though the idealized model differs from the physical ap¬ 
paratus, the orbit description within the idealized model 
is exact, and hence symplectic. For too coarse slicing, 
even though the orbit deviates from the “true” orbit, it 
remains “exact” (for a lattice with that slicing). One 
blames the inaccuracy on the fact that the lattice model 
deviates from the true lattice, not on the fact that the 
equations are being solved approximately. By varying 
the slicing one can judge from the limiting behavior what 
degree of slicing is appropriate. 

In a magnetic fields there is “geometric” horizontal fo¬ 
cusing even in a uniform field. (For example a cyclotron 
orbit with non-vanishing initial slope returns to its start¬ 
ing point after one full turn—its horizontal tune therefore 
being 1.0.) A deviant held index to introduces further 
horizontal and vertical focusing terms of equal magni¬ 
tude, but of opposite sign, into the focusing equations 

The corresponding focusing strength coefficients in 
electric fields are[l4] 


_ 

^x,l — 


o T 99 
^0 7o^o 
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TO 


(119) 


Note that for to = 0 (the “cylindrical” case) there is no 
vertical focusing. With the electric field being indepen¬ 
dent of y for cylindrical electrodes, this is obviously cor¬ 
rect. (The subscript “1” indicates “quadrupole” order; 


a more careful treatment brings in higher order correc¬ 
tions.) The parentheses segregate the geometric focusing 
(which is independent of to) from the gradient focusing 
(which vanishes for 1/r field variation). Re-writing these 
equations for our spherical to = 1 case; 


fc-sph. _ 
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'0 


2 2 
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( 120 ) 


With our sign convention, positive AT-value corresponds 
to focusing. Hence, for example, there is already “one 
unit” (in units of I/tq) of vertical focusing in the 
Coulomb case. In th e ETEAPOT formalism the focusing 
given by Eqs. ( fl^ is already implicitly included since 
the particle orbits are being tracked exactly. 

For modeling the focusing for to values other than to= 1 
the focusing coefficients we need are given by Eqs. (119). 


The geometric focusing terms already match but we have 
to make up the differences in gradient focusing: 


A _ rv(™) ivsph. _ m + l 
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( 121 ) 


These are the coefficients of the distributed focusing we 
have to apply to “convert” the Coulomb formulas to 
model the actual to value. The pure to= 1 Coulomb 
field itself provides “one unit” of vertical focusing. The 
strengths of the quads we insert artificially are therefore 
the sum of two terms; one cancels the Coulomb defo- 
cusing handled by the ETEAPOT formalism; the other 
superimposes focusing proportional to to, equal but op¬ 
posite horizontal and vertical. 

The artificial quadrupoles have zero length, but their 
length-strength product has to be matched to the “field 
integral” corresponding to length Abend of the bending 
slice being compensated. 

Before continuing with the treatment for to 7 ^ 0 it is 
important to remember the discontinuous increments to 
£ as a particle enters or leaves a bending element. The 
discontinuity is equal (in magnitude) to the change in 
potential energy. When correcting the focusing by artifi¬ 
cial quadrupoles we have to decide whether the artificial 
quadrupoles are “inside” or “outside”, since the actual 
deflections will be different in the two cases. Since the 
difference is quadratic in transverse displacement, the dif¬ 
ference would be called a “sextupole” effect. By treating 
the artificial quad as “inside”, which we do (because we 
decline to include longitudinal drift at zero potential), we 
avoid introducing an artificial sextupole effect. 

Whatever criteria for slicing quadrupoles that have 
been used for magnetic rings, the same criteria apply 
here. The default slicing in ETEAPOT treats a thick 
bend as a half bend, a single kick at the center, and then 
another half bend. For example, with to = —1.3, the cen¬ 
tral quad is vertically defocusing with inverse focal length 
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q = 1/f = 2.iilrQ, where £ is the arc length through the 
full bend. 

Numerically, with tq = 40 m and £ equal to, say, 
16 m, the ratio of focal length to element length is 
(40/16)^/2.3 ~ 2.7. With element length small compared 
to focal length, this suggests that the compensation will 
be fairly good even with such coarse slicing. The entire 
EDM ring, with its 8 full cells, would then be represented 
by 16 half-bends. In practice one will probably choose 
less coarse slicing. 

The transfer matrices for the thin effective quadrupole 
are 

=((-„ +Wo ?)' 

For m = 1 these kick matrices reduce to identity matri¬ 
ces. 

With linearized evolution through half bend for¬ 
mally represented by “transfer matrix” 
bend/kick/bend evolution through bend angle AO 
in an element with held index m is then described by 
the matrix 

(123) 

VIII. SPIN TRACKING 

r. Spin tracking through thick elements. The 

second main new feature implemented in ETEAPOT is 
spin evolution. To leading approximation spin precession 
occurs in central force, inverse square law, force held re¬ 
gions. With this assumption the orbit through a bend 
element of any individual particle lies in a single hxed 
plane. 

Each individual particle’s spin vector can be decom¬ 
posed into a (conserved) component s_l, normal to the 
bend plane, and a (processing) 2-component vector iy, 
lying in the bend plane. 

To lowest approximation hard edge bends are assumed, 
though a refractive dehection accompanying the change 
in electric potential is included. To a next approxima¬ 
tion the fringe held electric potential is represented as a 
linear ramp. In this region, because the particle speed 
is necessarily non-magic, there is a noticeable difference 
in precession rates of spin and momentum. Furthermore 
the precession error at entrance and exit add construc¬ 
tively. It is assumed that the paths through the entrance 
and exit fringe holds lie in the same plane as in the bend 
interior. 

As always in ETEAPOT, held deviations from radial 
inverse square law are modeled by artihcial quadrupole 
“kicks”. 

Quadrupoles, whether real or artihcial (as well as all 
other multipoles) are treated as thin. This includes the 


approximation that the electric potential is independent 
of position, both transversely and longitudinally in the 
element. Because the element thickness is treated as zero, 
there is no further time of hight error caused by the ne¬ 
glect of speed changes as the particle passes through the 
multipole. 

Superimposed multipoles are therefore exactly super¬ 
imposed which means that spin precession are modeled 
by successive rotations, with no dependence on the or¬ 
der in which they are applied. All such spin rotations 
are concatenated explicitly into a single (near-identity) 
precession matrix. For improved numerical precision all 
such concatenations are performed analytically and the 
result expressed as the exact sum of an identity matrix 
and a deviation matrix. This circumvents the problem 
of large, approximately canceling precessions and avoids 
spurious non-commutative geometric precessions. 

Quadrupoles too thick to be validly treated as thin, are 
sliced, with regions between the sliced thin quadrupoles 
treated as drifts. 

s. Spin representation. Laboratory frame spin 
components are [sx, Sy, Sz)- Bend frame spin components 
{sx,Sy,Sz), are illustrated on the left in Figure]^ which 
projects the spin vector onto the bend plane. 

fsx\ /-S||Sina\ 

“ j 1 ~ j ^-L 1 ■ (124) 

\S2/ \ S|| COSO J 

As with the orbit equation, the BMT equation is exactly 
solvable only because the orbit stays in a single bend 
plane. For particles in realistic beams this is very nearly, 
but not exactly, the local (x, y) central design plane. 

t. Transforming spin components from Lab 
frame to bend frame (and back again). From orbit 
tracking to some current location in the ring one knows 
the laboratory frame position, r, and momentum, p, and 
hence the angular momentum L = r x p. One also knows 
the spin vector s; 

r = Tx^i + TyJ + TzZ, 
p=PxX + pyy+pzZ, 

L = LxX -t- LyY + LzZ, 
s = Sxi +Syf + SzZ. (125) 

All of these quantities refer to a point infinitesimally past 
the bend entrance, after refractive correction, but not 
yet having included any entrance fringe held effect. Any 
spin precession in the interior (including fringe helds) of 
a bend element occurs in a single bend plane. To ex¬ 
ploit this reduction from 3D to 2D it is hrst necessary 
to obtain the spin components in an orthonormal frame 
having its “y” axis perpendicular to the plane and its 
“z” coordinate tangential to the orbit. The purpose of 
this section is to document this transformation. In the 
(always excellent) paraxial approximation, and for near¬ 
magic particle velocity, this transformation will be close 
to identity. 
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FIG. 3. (a) In the bend plane the spin vector s has precessed 
through angle a away from its nominal direction along the 
proton’s velocity. (Remember that different particles have dif¬ 
ferent bend planes.) (b) Projection of figure (a) onto the lab¬ 
oratory horizontal plane. The projected longitudinal axis is 
shown coinciding with the laboratory longitudinal axis, even 
if this is not exactly valid, x is the deviation of the (bold 
face) particle orbit from the (pale face) design orbit. If the 
bend plane coincides with the design bend plane (as is always 
approximately the case) Pq and z are identical. 9 is the ref¬ 
erence particle deviation angle from longitudinal and i9 is the 
tracked particle deviation angle from longitudinal. On the 
average 9 and i9 are the same, but betatron oscillations cause 
them to differ on a turn by turn basis, and also to make the 
instantaneous bend plane not quite horizontal. 


We can establish an orthonormal, right-handed basis 
triad with axis 3 parallel to p and axis 2 parallel to — L 
(where the negative sign is appropriate for clockwise or¬ 
bits); 


where R is an orthogonal matrix, 


/ an 

021 

031 

R = j ai 2 

022 

032 

\ai 3 

023 

033 


(130) 


(Aside: the magnitude |detR| of the determinant of R 
is necessarily 1, but the actual value is ±1. This sign 
correlates with the clockwise/counterclockwise orbit am¬ 
biguity.) 

Because R is orthogonal, R~^ = R^ and Eq. (129) 
can be inverted to give 


51 \ / 011 012 013 

5 2 1 = I 021 022 023 

,53/ \03i 032 O33 



(131) 


This yields the spin components in the bend frame. Their 
propagation through the bend is described below. 

After the bend plane spin components have been up¬ 
dated at the output of the bend element it is necessary 
to transform them back to the (local) laboratory frame. 
This entails repeating the preceeding formulas starting 
with Eqs (125), but with r, p and L having been eval¬ 


uated (in local laboratory coordinates) just inside the 
output face of the bend element. In other words all of 
the quantities in Eqs (125) will now refer to a point just 


before the bend exit. Then, following Eqs.(126) through 


(130) produces rotation matrix R which is now applicable 


at the bend output. Finally, modifying Eq. (131) appro¬ 


priately, the transformation back to laboratory compo¬ 
nents 


Px ^ ^ Py ^ ^ Pz ^ 
63 = —x-k—yH-z, 

P P P 

r X p 


- -L ’ 

ei = 02 X 63. 


(126) 


These equations can be re-expressed formally, with all 
coefficients known, as 


01 = oiix-k ai2y-k ai3Z 
02 = a2lX -k a22y + 0232 
03 = a3ix-k a32y-k a33Z. ( 127 ) 

The vector s can be expanded as 

S = Si 0 i -k 5262 + S303 

= si(aiix -k ai2y + ai3z) -k ... 

= (oiiSi + a2iS2 + a3iS3) X -k .... ( 128 ) 

The final relation can be expressed in matrix form as 



(129) 


( Sx\ fan 012 ai 3 \ /sA 

Sj, 1 = j 021 O22 O23 1 j S2 1 . (132) 

Sz / \ a 31 O32 O33 J yS3 J 

u. Exact solution of the BMT equation. We 

introduce, at least temporarily, the term “’’inverse square 
law bend” to characterize a bend having field index m = 
1, which is the case being treated in our 2D formalism. 
In this case the orbit stays in a single plane. Also, in 
this frame any precession of the spin is purely around an 
axis normal to the plane. Obtaining the initial values of 
the spin components in this frame was described in the 
previous section. 

In the bend plane the orbit lies in a single plane. Super¬ 
ficially this may suggest we are accounting only for hori¬ 
zontal betatron oscillations and neglecting vertical beta¬ 
tron oscillations. In fact, however, the ETEAPOT treat¬ 
ment accounts for arbitrary betatron and synchrotron 
motion by assigning different “wobbling planes” to each 
individual particle. Even allowing for vertical betatron 
motion these frames are all very nearly parallel to the 
global horizontal design frame of the ring. For the 2D 
evolution through electric bend elements in ETEAPOT, 
any betatron oscillations actually present for a particular 
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particle are folded into the determination of its particle- 
specific orbit plane, and the initial coordinates in this 
plane. 

As shown in Figure the initial spin vector is 

s = —S|| sinax-I-Sj^y-I-S|| cosaz. (133) 


from the known laboratory frame description of s, which 
also has to be updated as the particle exits the bend. 
(Ideally, in an EDM storage ring experiment any out-of- 
plane component of s would be evidence of non-vanishing 
electric dipole moment.) 

The precession rate of'd is governed by the equation 


Here s^y is the out-of-plane component of i, sy is the 
magnitude of the in-plane projection of s, and a is the 
angle between the projection of s onto the plane and the 
tangent vector to the orbit. 

Jackson’ s|l5) Eq. (11.171) gives the rate of change in 
an electric field E, of the longitudinal spin component as 


_d 

dt 


(/3-s) 


e 

rripC 


(s±,J • E) 




(134) 


(Note that Jackson’s sx.j is the component perpendic¬ 
ular to the tangent to the orbit not to the orbit plane.) 
Substituting from Eq. (133) the equation becomes 


-(»,cosa) = -— 


(135) 


With the orbit confined to a plane, any precession occurs 
about the normal to the plane, conserving Sy. Since the 
magnitude of s is conserved it follows that the magnitude 
S|| is also conse rved . This allo ws S| | to be treated as 
constant in Eq. (135). Then Eq. (135) reduces to 


da eE / g/3 1 \ 

dt rUpC \ 2 P J 


(136) 


dd d /s\ eE 
dt dt \r J p 


(137) 


where the curvature is 1/r = eE/(vp) and (just in this 
equation) s tempora rily s tands for arc length along the 
orbit. Dividing Eq. (136) by Eq. (137) and using pc = 


da 



(138) 


In this step we have also surrepticiously made the re¬ 
placement d —>■ 6. Though these angles are not the same, 
over arbitrarily long times they necessarily advance at the 
same rate. In any case the error in equating id to 9 be¬ 
comes progressively more valid in the fine-slicing limit, 
as the orbit is more nearly approximated by straight line 
segments. Explicitly the bend frame precession advance 
is the sum of two definite integrals 

Aa = — 1^ Ty — (139) 


where 


This is undoubtedly a fairly good approximate equation 
in any more-or-less constant electric field, but it is ex¬ 
act only for the m = 1 Keplerian electric field, which 
is the only field in which each orbit stays in its own 
fixed bend plane. In fact, the derivation is not quite 
valid even for our m = 1 case—though the design or¬ 
bit is circular, the betatron orbits are slightly elliptical. 
This violates our assumption that orbit and electric field 
are exactly orthogonal. Neglecting this error amounts to 
dropping a term from the RHS of Eq. (1341. This term 


is down by four orders of magnitude compared to the 
calculated value. Furthermore this term would average 
to zero (over many betatron cycles) except for a possible 
non-zero (quadratic in bend angle) commutation geomet¬ 
ric phase error. Such an error would be expected to be 
down by another four orders of magnitude. Furthermore, 
the importance of this “error” can be investigated, and 
reduced, by finer slicing of the bend element. 

Meanwhile the velocity vector itself has precessed by 
angle -d relative to a direction fixed in the laboratory. 
Note that this angle i), the angle of the individual parti¬ 
cle’s orbit is approximately, but not exactly equal to the 
angle 9 of the design orbit. 

In the ETEAPOT treatment each particle in an elec¬ 
tric bend element evolves in its own bend plane, sy is 
the component in this plane of the total spin vector. At 
every entrance to an electic bend iy has to be calculated 
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To account for fringe fields two more terms, Aa and 

-'FF,out 

Aa , w ill la ter be added directly to the right hand 
side of Eq. (139). This treatment makes the less-than- 


perfect assumption that the precession axes in the fringe 
fields are normal to the bend plane, which makes it valid 
to simply sum the bend field and fringe field precessions. 

V. Spin tracking through fringe fields. So far, 
as a particle enters or exits a bend element, its poten¬ 
tial energy has been treated as changing discontinuously 
with its kinetic energy changing correspondingly. For 
spin tracking, because there is significant excess spin pre¬ 
cession in fringe fields we have to treat this region more 
carefully. Instead of treating the potential as discontin¬ 
uous, we now assume the change occurs over a longi¬ 
tudinal distance Az^^ which, for estimation purposes, 
we take equal to the separation distance (symbol gap) 
between the electrodes; Az^^ = gap. (For the “proto- 
nium” model, introduced later for comparison with with 
analytic SCT calculations, the drift lengths are taken to 
be negligably small, making the fringe field spin preces¬ 
sion negligible.) For orbit tracking the fringe field region 
is assumed to be short enough to be treated as “thin”. 
That is, any change in the particle’s radial offset occuring 
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in range is to be neglected and the integrated de¬ 

flection applied at the center (i.e. the edge of the bend). 
As a result the curve x{z) is continuous, but its slope 
dx/dz is discontinuous. Entrance transitions from out¬ 
side a bend to inside are described first. 

Inside the bending element the increase in potential en¬ 
ergy from orbit centerline to radial position x is eAE [x). 
As synchroton oscillations move the particle radially in 
and out, the sign of AV{x), just inside the bend edge 
oscillates between negative and positive values, and the 
sign of the deviation from the magic velocity oscillates 
correspondingly. This will tend to average away the spin 
run-out occurring in the fringe field region over times long 
compared to the synchrotron period. In the long run it is 
the deviation from zero of this average that has to be de¬ 
termined. This can be done by pure numerical tracking 
or theoretically or, possibly, by a theoretically-weighted 
averaging of the numerical tracking data. 

Once one is able to determine the spin decoherence the 
task will shift to designing sextupole distributions capa¬ 
ble of increasing the spin coherence time SCT. Our ap¬ 
proach will be to study the effectiveness of such schemes 
before attempting to improve the precision of our fringe 
field treatment. 

The deflection angle of the design orbit in the 

fringe held at one such edge is approximately 

i 0.5 X 0.03/40 = 0.375 x lO'^); 

(141) 

this is half of the deflection occurring in advancing a dis¬ 
tance gap in the interior of the bend. (The angle is 

implicitly assumed to be positive, irrespective of whether 
the orbit is clockwise or counter-clockwise.) Consider a 
particle approaching the fringe held region at radial dis¬ 
placement X. At the longitudinal center of the fringe held 
region the kinetic energy of this particle deviates from its 
“proper” (i.e. fully-inside value at radial displacement x) 
by the amount 
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where AVtot is the total voltage increase from inner elec¬ 
trode to outer electrode. (The electric held points ra¬ 
dially inward in order for positive particles to bend to¬ 
ward negative x but, by convention, E is positive.) Here, 
for simplicity, we are neglecting the fact that the actual 
electric held will have more complicated x-dependence 
depending, for example, on the value of the held index 
m. Our assumed fringe held spatial dependence is also 
simplistic. 

The leading effect of passage through a bend region 
with 7 deviation from magic Ay, is a rate of change of 
spin angle a per unit dehection angle 0 given by 


Combining equations, the excess angular advance occur¬ 
ring while entering the bend at displacement x is 
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(Aside: it may be appropriate to keep another term in 
expansion (143) in order to include the fact that dis¬ 


persion introduces a correlation between 7 and x which, 
after averaging, leave s a hnite precession, even if < x > 
vanishes in Eq. ( 144[ ) .) 

In our initial treatment of this edge effect we are as¬ 
suming this precession lies in exactly the same plane as 

the orbit plane of the particle in the bend element, justi- 
-FF 

fying the notation Aa Entrance (and, later, exit) val¬ 
ues can then simply be added to the main precession 
through the bend element. Meanwhile, in the fringe field 
region the ad vance of the tangent to the orbit is as 

given by Eq. (141). The -I- sign on the rhs of Eq. (144) 


reflects the fact that, for a particle displaced radially out¬ 
ward, the particle momentum is completing some of its 
rotation in the fringe field where its magnitude is more 
positive than in the bend interior. 

Though the fringe field precession occurs continuously 
over the range gap it is applied discontinuously at the 
bend edge. This is consistent with our hard edge treat¬ 
ment of the particle’s momentum evolution. Because a is 
measured relative to the orbit direction, Eq. (144) gives 


the spin angle precession over and above the advance of 
the tangent to the orbit. 

The fact that spin and momentum angular advances 
do not match has come about because the particle has 
bent appreciably while its speed deviates from the magic 
value. On exiting the bend element the particle bends 
similarly while its 7 deviation is given by the same for¬ 
mula (142). Eq. (144) therefore applies to both entrance 


and exit. Unfortunately this means that excess input pre¬ 
cession and excess output precession combine construc¬ 
tively rather than tending to cancel (as edge focusing 
commonly does.) 

The largest magnitude can have is 

I A (FF)\ 1 E gap / e.g. 1 (10.5 x 10®) x 0.03 

’ ^ 4 m^cVe V 4 0.938 x 10® 


= 0.84 X 10"^.) 


(145) 


For a particle with magic velocity skimming the outer 
electrode, where the effect is maximum, the angular 
runout is given by 
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= 3.586 X (0.84 x 10”^) x (0.375 x 10"®) 


= 1.13 X 10 'radians/edge. 


(146) 
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+ 3.586 Ay.) (143) 


With perhaps 50 edges in the lattice, and revolution fre¬ 
quency of about 1 MHz, the maximum spin runout will 
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be about one revolution per second. This vastly exag¬ 
gerates the spin decoherence, of course, because it does 
not account for the averaging effect of synchrotron oscil¬ 
lations. A challenge for lattice design is to perfect the 
synchrotron oscillation averaging to zero. 

w. Spin tracking through thin elements. Most 
of the elements in a storage ring cause spin precession 
which approximately conserves the vertical component of 
spin Syf. The leading exceptions to this in a proton EDM 
storage rings are the quadrupoles (either focusing or de- 
focusing) present in the lattice to keep Py much smaller 
than fix. Particles having non-vanishing vertical beta¬ 
tron amplitude are deflected vertically which causes s^y 
to process. There is a very strong tendency for this pre¬ 
cession to cancel in subsequent quadrupoles and, there¬ 
fore, probably not contribute significantly to spin deco¬ 
herence. Nevertheless it is essential for this precession to 
be modeled correctly to avoid spurious EDM-mimicking 
precession. All quadrupoles and sextupoles in the lattice 
cause similar precession to at least some degree. 

In ETEAPOT the only thick elements are bends and 
drifts. Spin evolution through them has already been 
discussed. All other elements are treated as thin element 
(position dependent) kicks. Copying from the treatment 
in bends, spin evolution through a thin element is de¬ 
scribed by 


|Aa| a 



V2 7o / 


where the deflectionA0 is (for now) defined to be positive 
and Aa is the resulting angular deviation of the “bend” 
plane spin coordinate relative to the orbit; here “bend” 
is in quotes to emphasize the fact that the bend plane is 
not, in general, horizontal. By far the most important 
instance of this is the deflection of vertically offset orbits 
in lattice quadrupoles. For a particle with transverse 
position {x,y), the magnitude of the angular deviation 
A9 in a quadrupole of strength q is given by 

AO ^ \q\\/jp~+^, (148) 


where q is the inverse focal length of the quadrupole. The 
absolute value sign in this equation eventually has to be 
removed; it is included here so that the discussion of signs 
can be deferred. 

The magnitude of the angular deflections in sextupoles 
and higher order multipoles are a lso fu nctions only of the 
combination r = Eq. (1481 generalizes to 


A significant complication concerns the sign of Aa. For 
a particle with no vertical displacement there is no ambi¬ 
guity, since the bend plane in the quadrupole is the same 
as the overall (horizontal) lattice design plane. In this 
case, with y = 0, Eq. (147) can be made more explicit; 
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where, by conventio n, a h orizontal focusing quad has q > 
0. The sign in Eq. (150) reflects the fact that, for x > 
0 and g > 0, the quadrupole “helps” by bending the 
momentum in the same sense as the bending elements. 
This formula makes it clear that reversing the sign of q 
reverses the sign of Aa. 

For obtaining the proper sign for general multipoles 
it is necessary to handle consistently the transformation 
from laboratory to bend frame spin coordinates, which is 
why the sign issue has been deferred. 

For a particle incident at (a:o,yo) the equation of the 
line of intersection of the deflection plane with the trans¬ 
verse plane is 


y = yo- — iy- yo), 

Xo 


(151) 


In a quadrupole the roll-angle of the deflection plane 
(with counter-clockwise roll taken as positive) is (f>o = 
ta,n~^{yQ/xo), irrespective of quadrant and whether the 
quadrupole is focusing or defocusing. However the in¬ 
verse tangent function is, itself, multiple valued. To 
make it single valued one can determine (po using po = 

arctan2(g?/oj 9*0 )• 


Along with Eq. (147), this establishes both the sign 


and magnitude of Aa, while preserving the sign reversal 
when the sign of q reverses. Including also sextupoles, 
octupoles and other multipoles, one obtains 


<(’o,quad = larctan2(gyo,9a;o), 
<(>o.sext = 2arctan2(S'yo,S'a;o), 
4>o,oct = 3arctan2(Oi/o, Oa^o)- 


(152) 


With roll angles of the bend plane determined this way, 
the following formulas apply to all multipoles. 

The required spin rotation matrix is 

( Cll C12 Ci3 ' 

C21 C22 C23 
C3I C32 C33^ 

( cosAa 0 — sinAa 

0_ 1 0_ I RmiPo) 

sin Aa 0 cos Aa 


A6>quad = \q\r, A6>sext = A A6»oct = A (149) 

where S is the conventionally defined sextupole strength 
and O is the conventionally defined octupole strength. 


where 


( cosp 
sinp 
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— sin^ 0\ 

cost/) 0 

0 1 ) 


(153) 


^ The element strengths appearing in an SXF lattice description 


file include the integer factors. That is, the quad parameter 
is 61 = q, the sextupole parameter is 62 = S/2, the octupole 
parameter is 63 = S/6, and so on. 
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To avoid serious loss of numerical precision, it is essen¬ 
tial to use trigonometric identities to express the matrix 
product exactly and explicitly in this way. This utilizes 
the fact that the dominant, identity matrix part of the 
central matrix commutes with the outer matrices, whose 
product is I. The remaining, exact matrix elements are 


are explicitly small. 

Cii = — cos((/))^ sin(Aa/2), 

Ci 2 = cos{(j)) sin(Aa/2) sin((^), 

Ci 3 = — cos{(j)) cos(Aa;/2), 

C 21 = cos{cj)) sin(Aa/2) sin(0), 

C22 = -sin((/))^sin(Aa/2), 

C 23 = sin((/)) cos(Aa/2), 

C 31 = cos(0) cos(Aa/2), 

C32 = — sin(0) cos(Aq;/2), 

C33 = - sin(Aa/2) (154) 

Since the common factor 2sin(Aa/2) is tiny, all matrix 
elements are small even for angles (p of order tt. Multi¬ 
plying this matrix on the right by {sx, Sy, Sz)^ produces 
deviations {Asx, ^Sy, Asz) which, added to {sx, Sy, Sz), 
give the output spin coordinates. 
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